Goal:

Changed directory organization and adding git VC

Packages and user functions

knitr::opts_chunk$set(warning=F)
library(MCMCpack)
library(plotly) 
library(tidyverse)
library(driver) ## devtools::install_github("jsilve24/driver")
library(Rcpp)
library(ggsci)
library(cowplot)
#detach("package:vegan", unload=TRUE)

# functions -----------------------------------------------------------
# model prep
Sample_data <- function(Y, alpha, iter=1){
  d <- array(0, dim=c(N, D, iter))
  for (i in 1:nrow(Y)){
    d[i,,] <- t(rdirichlet(iter, Y[i,] + alpha[i,]))
  } 
  return(d)
}

# calculate quantiles
sourceCpp("~/Desktop/obmap/r_analysis/3dimOB/input/kzcolQuant.cpp")

#function for distance between 3d points
Dist3d <- function(df, row1, row2) {
  thing <- sqrt((df$AntPos[row1] - df$AntPos[row2])^2 +
                  (df$MedLat[row1] - df$MedLat[row2])^2 + 
                  (df$VenDor[row1] - df$VenDor[row2])^2)
  return(thing)
}

# 3d interactive scatterplot
Scat3d <- function(df, ml, ap, vd, color="voxrankperOR") {
  plot_ly(df, 
          x = ml, 
          y = ap, 
          z = vd, 
          color = color,
          text = ~paste('Gene:', olfrname, 
                        '<br>Rank:', voxrankperOR,
                        '<br>p50:', p50),
          marker = list(size = 6,
                        line = list(color = 'black',
                                    width = 0.5)),
          type = 'scatter3d',
          mode = 'markers') %>%
    layout(scene = list(xaxis = list(title = 'Medial-Lateral'),
                        yaxis = list(title = 'Anterior-Posterior'),
                        zaxis = list(title = 'Ventral-Dorsal')))
}

#given an OR, pick a number of high probability voxels based on signal to noise ratios
Scat_rank <- function(olfr, topX = 72, chooseOut = "plot", title = NA) {
  reranked <- ranked %>% 
    filter(olfrname == olfr) %>% 
    mutate(rankofrank = min_rank(voxRankSNR),
           rankcol = min_rank(desc(ifelse(rankofrank <= topX, rankofrank, NA))),
           isRanked = is.na(rankcol))
  
  besties <- reranked %>% filter(isRanked == 0) %>% arrange(voxRankSNR)
  worsties <- reranked %>% filter(isRanked == 1)
  all <- bind_rows(besties, worsties)
  
  if (chooseOut == "data") {
    return(all)
  } else {
    p <- plot_ly(type = "scatter3d", mode = "markers") %>% 
      add_trace(data=worsties, x=~AntPos, y=~MedLat, z=~VenDor, 
                color=~rankcol, opacity=0.15,
                text = ~paste('Gene:', olfrname, 
                              '<br>voxRankSNR:', voxRankSNR, 
                              '<br>voxSNRdim', voxSNRdim),
                marker = list(size = 6)) %>%
      add_trace(data=besties, x=~AntPos, y=~MedLat, z=~VenDor, color=~rankcol,
                text = ~paste('Gene:', olfrname, 
                              '<br>voxRankSNR:', voxRankSNR, 
                              '<br>voxSNRdim', voxSNRdim),
                marker = list(size = 6, line = list(color = 'black', width = 0.5))) %>%
      layout(title = title,
             scene = list(xaxis = list(title = 'Anterior-Posterior'),
                          yaxis = list(title = 'Medial-Lateral'),
                          zaxis = list(title = 'Ventral-Dorsal')))
    return(p)
  } #endif
} #endfunction

#plot p50 values for an olfrs top X voxels
p50plot <- function(olfr, rankx = 50, title = NA) {
  bestp50 <- ranked %>% 
    filter(olfrname == olfr) %>% 
    mutate(rankp1 = min_rank(p50), rankp = min_rank(desc(rankp1))) %>% filter(rankp <= rankx)
  worstp50 <- ranked %>% 
    filter(olfrname == olfr) %>% mutate(rankp = min_rank(p50)) %>% 
    filter(rankp > rankx) %>% mutate(rankna = NA)
  plot_ly(type = "scatter3d", mode = "markers") %>% 
    add_trace(data=worstp50, x=~AntPos, y=~MedLat, z=~VenDor, 
              color='rgb(10,10,10)', opacity=0.1,
              text = ~paste('Gene:', olfr, 
                            '<br>voxRankSNR:', voxRankSNR, 
                            '<br>voxSNRdim', voxSNRdim),
              marker = list(size = 6)) %>%
    add_trace(data=bestp50, x=~AntPos, y=~MedLat, z=~VenDor, color=~rankp,
              text = ~paste('Gene:', olfr, 
                            '<br>voxRankSNR:', voxRankSNR, 
                            '<br>voxSNRdim', voxSNRdim),
              marker = list(size = 6, line = list(color = 'black', width = 0.5))) %>%
    layout(title = title,
           scene = list(xaxis = list(title = 'Anterior-Posterior'),
                        yaxis = list(title = 'Medial-Lateral'),
                        zaxis = list(title = 'Ventral-Dorsal')))
} #endfunction

#define clusters of best X p50 points using a pairwise matrix and ability to output plots and data
#given a number of top ranking positions for an OR, cluster the points based on spatial position
Cluster <- function (olfr, topX = 72, minClustSize = 5, chooseOut = "plot", title = NA) {
  df <- ranked %>% 
    filter(olfrname == olfr) %>% 
    mutate(rankofrank = min_rank(desc(min_rank(p50)))) %>% 
    filter(rankofrank <= topX) %>% 
    arrange(desc(p50))
  
  cut <- ranked %>% 
    filter(olfrname == olfr) %>% 
    mutate(rankofrank = min_rank(desc(min_rank(p50)))) %>% 
    filter(rankofrank > topX) %>% 
    arrange(desc(p50))
  
  #make pairwise of matrixes
  dmatrix <- matrix(data = 0, nrow = nrow(df), ncol = nrow(df))
  for (i in 1:nrow(df)) {
    distances <- vector("numeric", length = nrow(df))
    toprankvox <- df$voxel[i]
    ap <- df$AntPos[i]
    ml <- df$MedLat[i]
    vd <- df$VenDor[i]
    for (j in 1:nrow(df)) {
      if (i == j) {
        distances[j] <- 0
      } else {
        distances[j] <- sqrt((df$AntPos[i] - df$AntPos[j])^2 + 
                               (df$MedLat[i] - df$MedLat[j])^2 +
                               (df$VenDor[i] - df$VenDor[j])^2)
      } #endif
    } #endforj
    neighbors <- which(distances <= sqrt(3))
    dmatrix[neighbors, i] <- 1
  } #endfori
  
  #cluster matrix
  cluster_list <- rep(NA, nrow(df))
  for (k in 1:ncol(dmatrix)) {
    #skip points that are already in a cluster
    if (is.na(cluster_list[k])) {
      round <- 1
      #do a bunch of rounds, find a way to have it run until it stops finding new points
      while (round < topX/3) {
        neigh <- which(dmatrix[,k] == 1)
        #for each neighbor of previous round of neighbors, find new neighbors
        for (l in 1:length(neigh)) {
          neigh <- c(neigh, which(dmatrix[,neigh[l]] == 1))
          neigh <- neigh[-which(duplicated(neigh))]
          neigh <- sort(neigh)
        } #endforl
        round <- round + 1
      } #endwhile
      cluster_list[neigh] <- k
    } else {
      next
    } #endif
  } #endfork
  
  df_out <- df %>% mutate(rawclust = cluster_list) %>% 
    group_by(rawclust) %>% 
    mutate(clustmaxp = max(p50), 
           clustminp = min(p50), 
           clustmeanp = mean(p50)) %>%
    add_tally() %>%
    ungroup() %>%
    mutate(clustmaxprank = dense_rank(desc(clustmaxp)), 
           clustmeanprank = dense_rank(desc(clustmeanp)), 
           clustsizerank = dense_rank(desc(n))) %>% 
    arrange(clustsizerank) %>%
    mutate(isCluster = ifelse(n >= minClustSize, 1, 0), 
           clust_unique = dense_rank(desc(clustmaxp * isCluster))) %>%
    select(p2.5:ORrankpervox, rankofrank, rawclust:clustsizerank, clust_unique, isCluster)
  
  df_clustered <- df_out %>% filter(isCluster == 1)
  too_small <- df_out %>% filter(isCluster == 0.5)
  
  cut_out <- cut %>% mutate(rawclust = NA,
                            clustmaxp = NA,
                            clustminp = NA,
                            clustmeanp = NA,
                            n = NA,
                            clustmeanprank = NA,
                            clustmaxprank = NA,
                            clustsizerank = NA,
                            isCluster = 0,
                            clust_unique = NA) %>%
    select(p2.5:ORrankpervox, 
           rankofrank, 
           rawclust:clustsizerank, 
           clust_unique, 
           isCluster) %>% 
    bind_rows(too_small)
  
  all_out <- bind_rows(df_clustered, cut_out) %>% unique()
  
  #return various things 'rgb(10,10,10)'
  if (chooseOut == "data") {
    return(all_out)
  } else {
    p <-  plot_ly(type = "scatter3d", mode = "markers") %>% 
      add_trace(data=cut_out, x=~AntPos, y=~MedLat, z=~VenDor, 
                color="shell", opacity=0.15,
                text = ~paste('Gene:', olfrname, 
                              '<br>C_size_rank:', clustsizerank, 
                              '<br>C_mean_p50:', clustmeanp,
                              '<br>C_max_p50:', clustmaxp,
                              '<br>Cluster:', clust_unique),
                marker = list(size = 6)) %>%
      add_trace(data=df_out, x=~AntPos, y=~MedLat, z=~VenDor, color=~clust_unique,
                text = ~paste('Gene:', olfrname, 
                              '<br>C_size_rank:', clustsizerank, 
                              '<br>C_mean_p50:', clustmeanp,
                              '<br>C_max_p50:', clustmaxp,
                              '<br>Cluster:', clust_unique),
                marker = list(size = 6, line = list(color = 'black', width = 0.5))) %>%
      layout(title = title,
             scene = list(xaxis = list(title = 'Anterior-Posterior'),
                          yaxis = list(title = 'Medial-Lateral'),
                          zaxis = list(title = 'Ventral-Dorsal')))
    return(p)
  } #endif
} #endfunction

#run cluster given an incrementing number of voxels to cluster from
#needs to also output some sort of summary statistic to define the optimal number of initial voxels that returns the "best" clusters 
BestML <- function(olfr, topMin = 50, topMax = 200, topBy = 25, minSize = 2,
                   clustersPerHalfBulb = 1, chooseOut = "plot", title = NA) {
  cphb <- 0
  topStep <- topMin
  while (cphb < clustersPerHalfBulb) {
    df_in <- Cluster(olfr, topX = topStep, minClustSize = minSize, chooseOut = "data")
    df_ml <- df_in %>% 
      filter(isCluster == 1) %>% 
      group_by(clust_unique) %>% 
      mutate(meanML = mean(MedLat),
             meanAP = mean(AntPos),
             meanVD = mean(VenDor),
             side = ifelse(meanML >= symline$mlvals[which(symline$apvals == round(meanAP))],
                           "Lateral", "Medial")) %>% 
      ungroup()
    
    df_bestM <- df_ml %>% 
      filter(side == "Medial") %>%
      mutate(MedRank = dense_rank(clust_unique),
             LatRank = NA,
             TopStep = topStep,
             sideRank = MedRank) %>%
      filter(MedRank <= clustersPerHalfBulb)
    
    df_bestL <- df_ml %>%
      filter(side == "Lateral") %>%
      mutate(MedRank = NA,
             LatRank = dense_rank(clust_unique),
             TopStep = topStep,
             sideRank = LatRank) %>%
      filter(LatRank <= clustersPerHalfBulb)
    
    clust_bestM <- unique(df_bestM$clust_unique)
    clust_bestL <- unique(df_bestL$clust_unique)
    
    #checks and counters
    cphb <- min(length(clust_bestM), length(clust_bestL))
    topStep <- topStep + topBy
  } #endwhile
  
  clust_bestML <- c(clust_bestM, clust_bestL)
  which_bestML <- df_in[-which(df_ml$clust_unique %in% clust_bestML),]
  
  df_notbest <- df_in %>% 
    mutate(sideRank = NA) %>% 
    select(-clust_unique) %>% 
    mutate(clust_unique = NA)
  df_best <- bind_rows(df_bestM, df_bestL)
  df_all <- bind_rows(df_best, df_notbest)
  
  #output
  if (chooseOut == "data") {
    return(df_all)
  } else if (chooseOut == "best") {
    return(df_best)
  } else if (chooseOut == "notbest") {
    return(df_notbest)
  } else {
    p <- plot_ly(type = "scatter3d", mode = "markers") %>% 
      add_trace(data=df_notbest, x=~AntPos, y=~MedLat, z=~VenDor, 
                color="shell", opacity=0.15,
                text = ~paste('Gene:', olfrname, 
                              '<br>C_size_rank:', clustsizerank, 
                              '<br>C_mean_p50:', clustmeanp,
                              '<br>C_max_p50:', clustmaxp,
                              '<br>Cluster:', clust_unique),
                marker = list(size = 6)) %>%
      add_trace(data=df_best, x=~AntPos, y=~MedLat, z=~VenDor, color=~clust_unique,
                text = ~paste('Gene:', olfrname, 
                              '<br>C_size_rank:', clustsizerank, 
                              '<br>C_mean_ML:', meanML,
                              '<br>C_max_p50:', clustmaxp,
                              '<br>Cluster:', clust_unique,
                              '<br>SideRank:', sideRank),
                marker = list(size = 6, line = list(color = 'black', width = 0.5))) %>%
      layout(title = title,
             scene = list(xaxis = list(title = 'Anterior-Posterior'),
                          yaxis = list(title = 'Medial-Lateral'),
                          zaxis = list(title = 'Ventral-Dorsal')))
    return(p)
  } #endif
} #endfunction

#given a list of Olfr names, output 1 medial and 1 lateral cluster for each name
#perhaps add an if or arg for really big lists to run a reduced step
ListML <- function(x, chooseOut = "plot", title = NA) {
  #use a list to build a df of unknown size instead of bind_row each iteration
  list_out <- vector("list", length = length(x))
  for (i in 1:length(x)) {
    list_out[[i]] <- BestML(x[i], topMin = 50, topMax = 150, topBy = 50, 
                            minSize = 2, clustersPerHalfBulb = 1, chooseOut = "best")
  } #endfor
  
  #always need notbest for the shape shell1
  notbest <- BestML(x[1], topMin = 100, topMax = 150, topBy = 50, minSize = 2, 
                    clustersPerHalfBulb = 1, chooseOut = "notbest")
  df_out <- bind_rows(list_out)
  
  #output
  if (chooseOut == "data") {
    return(df_out)
  } else if (chooseOut == "point") {
    df_point <- df_out %>% filter(p50 == clustmaxp)
    p <- plot_ly(type = "scatter3d", mode = "markers") %>% 
      add_trace(data=notbest, x=~AntPos, y=~MedLat, z=~VenDor, 
                color="shell", opacity=0.15,
                text = ~paste('Gene:', olfrname, 
                              '<br>C_size_rank:', clustsizerank, 
                              '<br>C_mean_p50:', clustmeanp,
                              '<br>C_max_p50:', clustmaxp,
                              '<br>Cluster:', clust_unique),
                marker = list(size = 6)) %>%
      add_trace(data=df_point, x=~AntPos, y=~MedLat, z=~VenDor, color=~olfrname,
                text = ~paste('Gene:', olfrname, 
                              '<br>C_size_rank:', clustsizerank, 
                              '<br>C_mean_ML:', meanML,
                              '<br>C_max_p50:', clustmaxp,
                              '<br>Cluster:', clust_unique,
                              '<br>SideRank:', sideRank),
                marker = list(size = 6, line = list(color = 'black', width = 0.5))) %>%
      layout(title = title,
             scene = list(xaxis = list(title = 'Anterior-Posterior'),
                          yaxis = list(title = 'Medial-Lateral'),
                          zaxis = list(title = 'Ventral-Dorsal')))
    return(p)
  } else {
    p <- plot_ly(type = "scatter3d", mode = "markers") %>% 
      add_trace(data=notbest, x=~AntPos, y=~MedLat, z=~VenDor, 
                color="shell", opacity=0.15,
                text = ~paste('Gene:', olfrname, 
                              '<br>C_size_rank:', clustsizerank, 
                              '<br>C_mean_p50:', clustmeanp,
                              '<br>C_max_p50:', clustmaxp,
                              '<br>Cluster:', clust_unique),
                marker = list(size = 6)) %>%
      add_trace(data=df_out, x=~AntPos, y=~MedLat, z=~VenDor, color=~olfrname,
                text = ~paste('Gene:', olfrname, 
                              '<br>C_size_rank:', clustsizerank, 
                              '<br>C_mean_ML:', meanML,
                              '<br>C_max_p50:', clustmaxp,
                              '<br>Cluster:', clust_unique,
                              '<br>SideRank:', sideRank),
                marker = list(size = 6, line = list(color = 'black', width = 0.5))) %>%
      layout(title = title,
             scene = list(xaxis = list(title = 'Anterior-Posterior'),
                          yaxis = list(title = 'Medial-Lateral'),
                          zaxis = list(title = 'Ventral-Dorsal')))
    return(p)
  } #endif
} #endfunction

#run Cluster and check if either top ranked M or L glom is dorsal and has known dorsal expression or is class 1. If True, find a dorsal glom for both M and L
DorsalML <- function(olfr, topMin = 100, topBy = 100, minSize = 2,
                     clustIn = 5, clustOut = 1, chooseOut = "plot") {
  print(olfr)
  clustFound <- 0
  topStep <- topMin
  while (clustFound != clustOut) {
    df_in <- Cluster(olfr, topX = topStep, minClustSize = minSize, chooseOut = "data")
    df_ml <- df_in %>% 
      filter(isCluster == 1) %>% 
      group_by(clust_unique) %>% 
      mutate(meanML = mean(MedLat),
             meanAP = mean(AntPos),
             meanVD = mean(VenDor)) %>%
      ungroup() %>%
      rowwise() %>%
      mutate(side = ifelse(meanML >= symline$mlvals[which(symline$apvals == round(meanAP))],
                           "Lateral", "Medial")) %>% 
      ungroup() %>%
      left_join(info, by = "olfrname") %>%
      rowwise() %>% 
      mutate(dorsalRating = sum(ifelse(class == 1, 2, 0), 
                                ifelse(tzsimple < 2, 1, 0),
                                ifelse(oe_region == "Dorsal", 1, 0),
                                na.rm = T)) %>%
      ungroup()
    
    df_bestM <- df_ml %>% 
      filter(side == "Medial") %>%
      mutate(MedRank = dense_rank(clust_unique),
             LatRank = NA,
             TopStep = topStep,
             minSize = minSize,
             clustIn = clustIn,
             sideRank = MedRank) %>%
      filter(MedRank <= clustIn) %>%
      arrange(MedRank)
    
    df_bestL <- df_ml %>%
      filter(side == "Lateral") %>%
      mutate(MedRank = NA,
             LatRank = dense_rank(clust_unique),
             TopStep = topStep,
             minSize = minSize,
             clustIn = clustIn,
             sideRank = LatRank) %>%
      filter(LatRank <= clustIn) %>%
      arrange(LatRank)
    
    max <- max(c(df_bestM$meanVD[1], df_bestL$meanVD[1]))
    min <- min(c(df_bestM$meanVD[1], df_bestL$meanVD[1]))
    
    #would be nice to remake using dplyr::case_when()
    if (is.na(max)) {
      print(paste("NA", as.character(topStep)))
      topStep <- topStep + topBy
      next
    } else {
      if (max >= 13) {
        if (min < 13) {
          if (df_bestM$dorsalRating[1] >= 2) {
            #OR has dorsal info, constrain both med/lat to have dorsal glom
            if (df_bestM$meanVD[1] == max) {
              #Medial was dorsal, find dorsal lateral glom and viceversa
              df_bestM <- df_bestM %>%
                filter(MedRank == 1) %>% mutate(test = "1m")
              df_bestL <- df_bestL %>%
                filter(meanVD >= 13) %>%
                filter(LatRank == min(LatRank)) %>% mutate(test = "1m")
              
            } else if (df_bestL$meanVD[1] == max) {
              df_bestL <- df_bestL %>%
                filter(LatRank == 1) %>% mutate(test = "1l")
              df_bestM <- df_bestM %>% 
                filter(meanVD >= 13) %>%
                filter(MedRank == min(MedRank)) %>% mutate(test = "1l")
            } #endif df_bestM$meanVD[1] == max
          } else if (df_bestM$dorsalRating[1] == 1) {
            #OR has mixed info, return best
            df_bestM <- df_bestM %>% filter(MedRank == 1)  %>% mutate(test = "2")
            df_bestL <- df_bestL %>% filter(LatRank == 1)  %>% mutate(test = "2")
          } else {
            #OR has ventral info, constrain both med/lat to have ventral glom
            if (df_bestM$meanVD[1] == min) {
              #Medial was ventral, find ventral ateral glom and viceversa
              df_bestM <- df_bestM %>%
                filter(MedRank == 1)   %>% mutate(test = "3m")
              df_bestL <- df_bestL %>%
                filter(meanVD < 13) %>%
                filter(LatRank == min(LatRank))  %>% mutate(test = "3m")
              
            } else if (df_bestL$meanVD[1] == min) {
              df_bestL <- df_bestL %>%
                filter(LatRank == 1)  %>% mutate(test = "3l")
              df_bestM <- df_bestM %>% 
                filter(meanVD < 13) %>%
                filter(MedRank == min(MedRank))  %>% mutate(test = "3l")
            } #endif (df_bestM$meanVD[1] == min)
          } #endif (df_bestM$dorsalRating[1] >= 2)
        } #endif (min < 13)
      } else {
        #both gloms are ventral, check if dorsal info
        if (df_bestM$dorsalRating[1] >= 2) {
          #both gloms are ventral, but OR has dorsal info, constrain to dorsal
          df_bestM <- df_bestM %>%
            filter(meanVD >= 13) %>%
            filter(MedRank == min(MedRank))  %>% mutate(test = "4")
          df_bestL <- df_bestL %>%
            filter(meanVD >= 13) %>%
            filter(LatRank == min(LatRank))  %>% mutate(test = "4")
        } else {
          df_bestM <- df_bestM %>% filter(MedRank == 1)   %>% mutate(test = "5")
          df_bestL <- df_bestL %>% filter(LatRank == 1)  %>% mutate(test = "5")
        } #endif (df_bestM$dorsalRating[1] >= 2)
      } #endif (max >= 13)
    } #endif (is.nax(max))
    
    df_bestM <- df_bestM %>% filter(MedRank == min(MedRank))  %>% mutate(test = "suck")
    df_bestL <- df_bestL %>% filter(LatRank == min(LatRank)) %>% mutate(test = "suck")
    
    clust_bestM <- unique(df_bestM$clust_unique)
    clust_bestL <- unique(df_bestL$clust_unique)
    
    #checks and counters
    clustFound <- min(c(length(clust_bestM), length(clust_bestL)))
    print(topStep)
    topStep <- topStep + topBy
    clustIn <- round(clustIn * 1.5)
  } #endwhile
  
  clust_bestML <- c(clust_bestM, clust_bestL)
  which_bestML <- df_in[-which(df_ml$clust_unique %in% clust_bestML),]
  
  df_best <- bind_rows(df_bestM, df_bestL) %>% unique()
  df_notbest <- df_in %>% 
    mutate(sideRank = NA) %>% 
    select(-clust_unique) %>% 
    mutate(clust_unique = NA) %>%
    unique()
  df_all <- bind_rows(df_best, df_notbest) %>% unique()
  
  #output
  if (chooseOut == "data") {
    return(df_all)
  } else if (chooseOut == "best") {
    return(df_best)
  } else if (chooseOut == "notbest") {
    return(df_notbest)
  } else {
    p <- plot_ly(type = "scatter3d", mode = "markers") %>%
      add_trace(data=df_notbest, x=~AntPos, y=~MedLat, z=~VenDor,
                color="shell", opacity=0.15,
                text = ~paste('Gene:', olfrname,
                              '<br>C_size_rank:', clustsizerank,
                              '<br>C_mean_p50:', clustmeanp,
                              '<br>C_max_p50:', clustmaxp,
                              '<br>Cluster:', clust_unique),
                marker = list(size = 6)) %>%
      add_trace(data=df_best, x=~AntPos, y=~MedLat, z=~VenDor, color=~clust_unique,
                text = ~paste('Gene:', olfrname,
                              '<br>C_size_rank:', clustsizerank,
                              '<br>C_mean_ML:', meanML,
                              '<br>C_max_p50:', clustmaxp,
                              '<br>Cluster:', clust_unique,
                              '<br>SideRank:', sideRank),
                marker = list(size = 6, line = list(color = 'black', width = 0.5))) %>%
      layout(title = title,
             scene = list(xaxis = list(title = 'Anterior-Posterior'),
                          yaxis = list(title = 'Medial-Lateral'),
                          zaxis = list(title = 'Ventral-Dorsal')))
    return(p)
  } #endif
} #endfunction

#given a list of Olfr names, output 1 medial and 1 lateral cluster for each name
ListDorML <- function(x, chooseOut = "plot", title = NA) {
  #use a list to build a df of unknown size instead of bind_row each iteration
  list_out <- vector("list", length = length(x))
  for (i in 1:length(x)) {
    list_out[[i]] <- DorsalML(x[i], topBy = 200, chooseOut = "best")
  } #endfor
  
  #always need notbest for the shape shell1
  notbest <- DorsalML(x[1], chooseOut = "notbest")
  df_out <- bind_rows(list_out)
  
  #output
  if (chooseOut == "data") {
    return(df_out)
  } else if (chooseOut == "point") {
    df_point <- df_out %>% filter(p50 == clustmaxp)
    p <- plot_ly(type = "scatter3d", mode = "markers") %>% 
      add_trace(data=notbest, x=~AntPos, y=~MedLat, z=~VenDor, 
                color="shell", opacity=0.15,
                text = ~paste('Gene:', olfrname, 
                              '<br>C_size_rank:', clustsizerank, 
                              '<br>C_mean_p50:', clustmeanp,
                              '<br>C_max_p50:', clustmaxp,
                              '<br>Cluster:', clust_unique),
                marker = list(size = 6)) %>%
      add_trace(data=df_point, x=~AntPos, y=~MedLat, z=~VenDor, color=~olfrname,
                text = ~paste('Gene:', olfrname, 
                              '<br>C_size_rank:', clustsizerank, 
                              '<br>C_mean_ML:', meanML,
                              '<br>C_max_p50:', clustmaxp,
                              '<br>Cluster:', clust_unique,
                              '<br>SideRank:', sideRank),
                marker = list(size = 6, line = list(color = 'black', width = 0.5))) %>%
      layout(title = title,
             scene = list(xaxis = list(title = 'Anterior-Posterior'),
                          yaxis = list(title = 'Medial-Lateral'),
                          zaxis = list(title = 'Ventral-Dorsal')))
    return(p)
  } else {
    p <- plot_ly(type = "scatter3d", mode = "markers") %>% 
      add_trace(data=notbest, x=~AntPos, y=~MedLat, z=~VenDor, 
                color="shell", opacity=0.15,
                text = ~paste('Gene:', olfrname, 
                              '<br>C_size_rank:', clustsizerank, 
                              '<br>C_mean_p50:', clustmeanp,
                              '<br>C_max_p50:', clustmaxp,
                              '<br>Cluster:', clust_unique),
                marker = list(size = 6)) %>%
      add_trace(data=df_out, x=~AntPos, y=~MedLat, z=~VenDor, color=~olfrname,
                text = ~paste('Gene:', olfrname, 
                              '<br>C_size_rank:', clustsizerank, 
                              '<br>C_mean_ML:', meanML,
                              '<br>C_max_p50:', clustmaxp,
                              '<br>Cluster:', clust_unique,
                              '<br>SideRank:', sideRank),
                marker = list(size = 6, line = list(color = 'black', width = 0.5))) %>%
      layout(title = title,
             scene = list(xaxis = list(title = 'Anterior-Posterior'),
                          yaxis = list(title = 'Medial-Lateral'),
                          zaxis = list(title = 'Ventral-Dorsal')))
    return(p)
  } #endif
} #endfunction

#using heatmap peak data, create a 3d raw point or find nearest blank OB shell point to that raw point, options to use tan zone to assign a VD position or use average of VD peaks
Heat3D <- function(olfr, dimrep = 1, dv = "oe", raw = F, chooseOut = "plot") {
  df_in <- heatmap_peaks %>% filter(olfrname %in% olfr) %>% arrange(side)
  
  if (nrow(df_in) == 0) {
    print(paste("No", olfr, "Found"))
    next
  } #endif
  
  if (dv == "oe") {
    df_dv <- df_in %>%
      mutate(VD_selected = VDtz,
             VD_value = "Tan Based")
  } else {
    df_dv <- df_in %>%
      mutate(VD_selected = VDavg,
             VD_value = "Avg Based")
  } #endif
  
  if (dimrep == 1) {
    df_out <- df_dv %>% 
      select(olfrname, AntPos10, APval10, MedLat8, MLval8, VD_selected, VD_value, side) %>%
      rename(AP_selected = AntPos10, AP_value = APval10,
             ML_selected = MedLat8, ML_value = MLval8)
  } else if (dimrep == 2) {
    df_out <- df_dv %>% 
      select(olfrname, AntPos13, APval13, MedLat11, MLval11, VD_selected, VD_value, side) %>%
      rename(AP_selected = AntPos13, AP_value = APval13,
             ML_selected = MedLat11, ML_value = MLval11)
  } else if (dimrep == 3) {
    df_out <- df_dv %>% 
      select(olfrname, AntPos15, APval15, MedLat16, MLval16, VD_selected, VD_value, side) %>%
      rename(AP_selected = AntPos15, AP_value = APval15,
             ML_selected = MedLat16, ML_value = MLval16)
  } else if (dimrep == "123") {
    df_out <- df_dv %>% rowwise() %>%
      mutate(AP_selected = mean(c(AntPos10, AntPos13, AntPos15)),
             AP_value = mean(c(APval10, APval13, APval15)),
             ML_selected = mean(c(MedLat8, MedLat11, MedLat16)),
             ML_value = mean(c(MLval8, MLval11, MLval16))) %>%
      ungroup()
  } else if (dimrep == "12") {
    df_out <- df_dv %>% rowwise() %>%
      mutate(AP_selected = mean(c(AntPos10, AntPos13)),
             AP_value = mean(c(APval10, APval13)),
             ML_selected = mean(c(MedLat8, MedLat11)),
             ML_value = mean(c(MLval8, MLval11))) %>%
      ungroup()
  } #endif
  
  if (raw == T) {
    if (chooseOut == "plot") {
      plot <- plot_ly(type = "scatter3d", mode = "markers") %>%
        add_trace(data = blankdata,
                  x = ~AntPos, y=~MedLat, z=~VenDor,
                  color="shell", opacity=0.15,
                  text = ~paste('AntPos:', AntPos,
                                '<br>MedLat:', MedLat,
                                '<br>VenDor:', VenDor),
                  marker = list(size = 6)) %>%
        add_trace(data = df_out, 
                  x = ~AP_selected, 
                  y = ~ML_selected, 
                  z = ~VD_selected, 
                  color = ~olfrname,
                  text = ~paste('Gene:', olfrname, 
                                '<br>AntPos:', AP_selected,
                                '<br>APvalue', AP_value,
                                '<br>MedLat:', ML_selected,
                                '<br>MLvalue:', ML_value,
                                '<br>VenDor:', VD_selected,
                                '<br>VDvalue:', VD_value,
                                '<br>Side:', side,
                                '<br>Point:', "Raw"),
                  marker = list(size = 6, 
                                line = list(color='black', width = 0.5))) %>%
        layout(scene = list(xaxis = list(title = 'Anterior-Posterior'),
                            yaxis = list(title = 'Medial-Lateral'),
                            zaxis = list(title = 'Ventral-Dorsal')))
      return(plot)
    } else if (chooseOut == "data") {
      return(df_out)
    }
  } else if (raw == F) { 
    df_close <- tibble("AntPos" = numeric(), "MedLat" = numeric(), "VenDor" = numeric())
    for (i in 1:nrow(df_out)) {
      distances <- vector(mode = "numeric", length = nrow(blankdata))
      sidevec <- c("Lateral", "Medial")
      for (j in 1:nrow(blankdata)) {
        distances[j] <- sqrt((df_out$AP_selected[i] - blankdata$AntPos[j])^2 +
                               (df_out$ML_selected[i] - blankdata$MedLat[j])^2 + 
                               (df_out$VD_selected[i] - blankdata$VenDor[j])^2)
        df_close[i,] <- blankdata[which(rank(distances, ties.method = "first") == 1),]
      }
    }
    
    df_close_out <- df_close %>% 
      mutate(side = sidevec) %>%
      left_join(df_out, by = "side") %>%
      select(olfrname, AntPos, VenDor, MedLat, side, everything())
    
    if (chooseOut == "plot") {
      plot <- plot_ly(type = "scatter3d", mode = "markers") %>%
        add_trace(data = blankdata,
                  x = ~AntPos, y=~MedLat, z=~VenDor,
                  color="shell", opacity=0.15,
                  text = ~paste('AntPos:', AntPos,
                                '<br>MedLat:', MedLat,
                                '<br>VenDor:', VenDor),
                  marker = list(size = 6)) %>%
        add_trace(data = df_close_out, 
                  x = ~AntPos, 
                  y = ~MedLat, 
                  z = ~VenDor, 
                  color = ~olfrname,
                  text = ~paste('Gene:', olfrname, 
                                '<br>AntPos:', AP_selected,
                                '<br>APvalue', AP_value,
                                '<br>MedLat:', ML_selected,
                                '<br>MLvalue:', ML_value,
                                '<br>VenDor:', VD_selected,
                                '<br>VDvalue:', VD_value,
                                '<br>Side:', side,
                                '<br>Point:', "Nearest Shell"),
                  marker = list(size = 6, 
                                line = list(color='black', width = 0.5))) %>%
        layout(scene = list(xaxis = list(title = 'Anterior-Posterior'),
                            yaxis = list(title = 'Medial-Lateral'),
                            zaxis = list(title = 'Ventral-Dorsal')))
      return(plot)
    } else if (chooseOut == "data") {
      return(df_close_out)
    } #endif
  } #endif
} #endfunction

#compute distances between heatmap peak point and good point data featuring eurovision movie memes
DistHeat3D <- function(olfr_list, heat_dimrep = 1, heat_dv = "oe", heat_raw = F) {
  doubletrouble <- tibble(olfrname = character(), side = character(), 
                          distance = numeric(), ap_ml_vd_3d = character(), 
                          ap_ml_vd_heat = character(), .rows = 0)
  
  for (i in 1:length(olfr_list)) {
    heatpoint <- Heat3D(olfr_list[i], dimrep = heat_dimrep, 
                        dv = heat_dv, raw = heat_raw, chooseOut = "data") %>%
      select(olfrname, AntPos, MedLat, VenDor, side) %>% 
      mutate(origin = "heat",
             ap_ml_vd_heat = paste(AntPos, MedLat, VenDor,sep = "_"))
    
    goodpoint <- good_point %>% 
      filter(olfrname == olfr_list[i]) %>% 
      select(olfrname, AntPos, MedLat, VenDor, side) %>% 
      mutate(origin = "3d",
             ap_ml_vd_3d = paste(AntPos, MedLat, VenDor,sep = "_")) 
    
    doubletrouble <- bind_rows(doubletrouble, 
                                  left_join(heatpoint, goodpoint, 
                                            by = c("olfrname", "side"), 
                                            suffix = c("_h", "_3")) %>%
                                    mutate(distance = sqrt((AntPos_h - AntPos_3)^2 +
                                                             (MedLat_h - MedLat_3)^2 +
                                                             (VenDor_h - VenDor_3)^2)) %>%
                                    select(olfrname, side, distance, 
                                           ap_ml_vd_3d, ap_ml_vd_heat))
  } #endfor
  doubletrouble <- doubletrouble %>% 
    mutate(dimrep = heat_dimrep,
           vd_assignment = heat_dv,
           use_raw = heat_raw)
  return(doubletrouble)
} #endfunction

Build model

Input: TPM (gene length and sequence depth normalized expression unit) for 1088 intact ORs across 269 samples across 6 Anterior-Posterior, 3 Medial-Lteral, and 3 Ventral-Dorsal 100um vibratome sliced single OBs.

kzY <- read_csv("~/Desktop/obmap/r_analysis/3dimOB/input/covarintactchemo_over_samples_200923.csv", col_names = TRUE) %>% select(-X1) %>% filter(dimrep == 1 | dimrep == 2)
ectopic <- c("Olfr287","Olfr32")
kzYgood <- dplyr::select(kzY, -ectopic) %>% dplyr::select(-name, -rep, -slice, -dim, -dimrep)
dim(kzYgood) #1100 minus number of genes in ectopic
kzmY <- as.matrix(kzYgood)
rownames(kzmY) <- paste0("sample", 1:dim(kzY)[1])

#rdirichlet requires vector not dataframe so matrix the tibble
kzX <- read_csv("~/Desktop/obmap/r_analysis/3dimOB/input/allmice_covariates_trim_voxweights_v3.csv") %>% filter(name %in% kzY$name)
N <- dim(kzmY)[1] #number of sections
D <-  dim(kzmY)[2] #number of genes

# Voxel System
d1 <- 23 #antpost
d2 <- 22 #medlat
d3 <- 23 #vendor

#load voxels from bulb_in_cube 
# voxalls_test <- read_csv("~/Desktop/obmap/r_analysis/3dimOB/input/blankOBcoords200922.csv") %>% mutate(type = "latfill") #incorporating fill for 1 voxel thick layer on most lateral surface
# voxalls_0 <- readRDS("~/Desktop/obmap/r_analysis/3dimOB/input/191010_outer_inner_coords_ml23nohole.RDS")
# voxalls_0 <- bind_rows(voxalls_0, voxalls_test)
# voxalls_1 <- voxalls_0[-which(duplicated(voxalls_0)),]
# vox_ap <- voxalls_1[which(voxalls_1$AntPos <= d1),]
# vox_ml <- vox_ap[which(vox_ap$MedLat <= d2),]
# vox_vd <- vox_ml[which(vox_ml$VenDor <= d3),]
# voxalls <- vox_vd %>% select(-type) %>% unique()
voxalls <- read_csv("~/Desktop/obmap/r_analysis/3dimOB/input/voxalls_200924.csv")

# voxweights <- voxalls %>% 
#   group_by(AntPos) %>% 
#   mutate(ap.cnt = n(), ap.weight = ap.cnt/dim(voxalls)[1]) %>% 
#   ungroup() %>%
#   group_by(VenDor) %>% 
#   mutate(vd.cnt = n(), vd.weight = vd.cnt/dim(voxalls)[1]) %>% 
#   ungroup() %>%
#   group_by(MedLat) %>% 
#   mutate(ml.cnt = n(), ml.weight = ml.cnt/dim(voxalls)[1]) %>% 
#   ungroup() %>%
#   mutate(voxweight = ap.weight*vd.weight*ml.weight)
# 
# vw_ap <- voxalls %>% 
#   group_by(AntPos) %>% 
#   summarise(ap.cnt = n(), ap.weight = ap.cnt/dim(voxalls)[1], dim = "AntPos") %>% 
#   rename(slice = AntPos, count = ap.cnt, weight = ap.weight)
# vw_ml <- voxalls %>% 
#   group_by(MedLat) %>% 
#   summarise(ml.cnt = n(), ml.weight = ml.cnt/dim(voxalls)[1], dim = "MedLat") %>% 
#   rename(slice = MedLat, count = ml.cnt, weight = ml.weight)
# vw_vd <- voxalls %>% 
#   group_by(VenDor) %>% 
#   summarise(vd.cnt = n(), vd.weight = vd.cnt/dim(voxalls)[1], dim = "VenDor") %>% 
#   rename(slice = VenDor, count = vd.cnt, weight = vd.weight)
# voxweights <- bind_rows(vw_ml, vw_vd, vw_ap)
#write_csv(voxweights, "~/Desktop/obmap/r_analysis/3dimOB/input/voxweights_200923.csv")
voxweights <- read_csv("~/Desktop/obmap/r_analysis/3dimOB/input/voxweights_200923.csv")

# pick prior and Monte Carlo Sample 
alpha <- matrix(0.65, N, D)
iter <- 200
pi_post <- Sample_data(kzmY, alpha, iter) 

check <- ifelse(t(voxalls[,kzX$dim]) == kzX$slice, T, F)
weights <- ifelse(check, kzX$weight, 0)
weights <- miniclo(t(weights)) 
colnames(weights) <- rownames(kzmY) 

# We are going to use the ILR basis for computation
contrast.matrix <- create_default_ilr_base(D)
eta_post <- ilr_array(pi_post, contrast.matrix, 2)

# now calculate weighted composition
w8_start <- Sys.time()
q_ilr <- array(0, dim=c(nrow(voxalls), D-1, iter))
for (i in 1:iter) q_ilr[,,i] <- weights %*% eta_post[,,i] #matrix multiplication operator
w8_end <- Sys.time()
(w8_time <- w8_end - w8_start)

#clear RAM
remove(pi_post)
remove(eta_post)
gc()

#proportional composition rather than in ILR, using driver 18-08-16 new ilrInv_array with transposed q_ilr and coord == 1
propstart <- Sys.time()
ilr_prop <- aperm(q_ilr, c(2,1,3))

remove(q_ilr)
gc() 
q_proportions <- ilrInv_array(ilr_prop, contrast.matrix, 1)
propend <- Sys.time()
(proptime <- propend - propstart)

remove(ilr_prop)
remove(contrast.matrix)
gc()

sumstr <- Sys.time()
q_proptidy <- aperm(q_proportions, c(3,2,1))
d <- dim(q_proptidy)
q_proptidy <- matrix(q_proptidy, d[1], prod(d[-1]))
quantiles <- c(0.025, .25, .5, .75, .975)
q_prop_quantiles <- t(kzcolQuant(q_proptidy, quantiles))
colnames(q_prop_quantiles) <- paste0("p", 100*quantiles)
voxel=rep(1:d[2], times=D)
coord.ap = voxalls[voxel, "AntPos"] 
coord.ml = voxalls[voxel, "MedLat"]
coord.vd = voxalls[voxel, "VenDor"]
tidy_result <- data.frame(voxel=rep(1:d[2], times=D), 
                          gene=rep(1:d[3], each=nrow(voxalls)), 
                          q_prop_quantiles) %>% 
  bind_cols(coord.ap, coord.ml, coord.vd)
head(tidy_result)
sumend <- Sys.time()
(sumtime <- sumend - sumstr)
saveRDS(tidy_result, file = "~/Desktop/obmap/r_analysis/3dimOB/tidyresults/200i_allchemo_dimrep12_allvox.RDS")

Load info and features

Both sexes represented in each dimensional group. Line of glomeruli symmetry determined from 1D AP,ML heatmap mean of top 2 peaks Olfr features - class, tanzone, matsunami OE diffE, RTPdko OB scaffold with 2922 voxels, OR positions will only be assigned to these positions

kzY <- read_csv("~/Desktop/obmap/r_analysis/3dimOB/input/covarintactchemo_over_samples_200923.csv", col_names = TRUE) %>% select(-X1) %>% filter(dimrep == 1 | dimrep == 2)
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   name = col_character(),
##   dim = col_character()
## )
## See spec(...) for full column specifications.
ectopic <- c("Olfr287","Olfr32")
kzYgood <- dplyr::select(kzY, -ectopic) %>% dplyr::select(-name, -rep, -slice, -dim, -dimrep)
## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(ectopic)` instead of `ectopic` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
dim(kzYgood) #1100 minus number of genes in ectopic
## [1]  135 1482
kzmY <- as.matrix(kzYgood)
rownames(kzmY) <- paste0("sample", 1:dim(kzY)[1])

#dimwise TPM plot
s1 <- kzY %>% select(slice, dim, dimrep, Olfr1377) %>% mutate(olfr = "1377") %>% rename("gene" = "Olfr1377")
s8 <- kzY %>% select(slice, dim, dimrep, Olfr881) %>% mutate(olfr = "881") %>% rename("gene" = "Olfr881")
scombo <- bind_rows(s1,s8)
ggplot(scombo) + 
  geom_line(aes(slice, gene, color = olfr)) + 
  facet_grid(vars(dimrep), vars(dim), scales = "free_y") +
  ggtitle("TPM of Olfr881/1377 across dimensions and replicates")

#Examine Olfr881 TPM for all mice replicates from all dims
ggplot(kzY) + geom_line(aes(slice, Olfr881)) + 
  facet_grid(vars(dimrep), vars(dim), scales = "free_y") +
  ggtitle("TPM of Olfr881 across dimensions and replicates")

#feature info
info <- read_csv("~/Desktop/obmap/r_analysis/3dimOB/input/knowntanwavgFI.csv", col_names = TRUE) %>% 
  rename("olfrname" = "gene") %>%
  select(olfrname:RTP, known, lowTPM)
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   gene = col_character(),
##   tan_zone = col_character(),
##   oe_region = col_character(),
##   RTP = col_character(),
##   fisurface = col_logical(),
##   ish_id = col_logical(),
##   lacz_xg = col_logical(),
##   fp_xg = col_logical(),
##   med_glom = col_character(),
##   lat_glom = col_character(),
##   citation = col_character(),
##   known = col_logical(),
##   lowTPM = col_logical()
## )
## See spec(...) for full column specifications.
#line of bulb symmetry as found using single-dimension heatmap data
#the average of the two calculated lines is used to call whether a predicted glomeruli position is medial or lateral
symline <- read_csv("~/Desktop/obmap/r_analysis/3dimOB/input/symline.csv")
## Parsed with column specification:
## cols(
##   apvals = col_double(),
##   mlvals = col_double()
## )
ggplot() +
  geom_blank() +
  geom_abline(aes(slope = 0.32655, intercept = 7.61119, color = "rep1")) +
  geom_abline(aes(slope = 0.24922, intercept = 8.83890, color = "rep2")) +
  geom_abline(aes(slope = 0.287885, intercept = 8.225045, color = "mean")) +
  xlim(0,23) +
  ylim(0,22) + 
  ggtitle("Glomeruli symmetry line as calculated by single dimension heatmap peaks") +
  xlab("Mean position of top 2 AP peaks") +
  ylab("Mean position of top 2 ML peaks")

# Example 3D left OB for future orientation of dimensions
#blankdata <- Scat_rank("Olfr881", 1, "data") %>% select(AntPos:VenDor) %>% arrange(AntPos, MedLat, VenDor)
#write_csv(blankdata, "~/Desktop/rproj/obmap/allmice/v21_gen25/base_files/blankOBcoords.csv")
#blankOBcoords200922
blankdata <- read_csv("~/Desktop/obmap/r_analysis/3dimOB/input/blankOBcoords200922.csv")
## Parsed with column specification:
## cols(
##   AntPos = col_double(),
##   MedLat = col_double(),
##   VenDor = col_double()
## )
dvml_blank <- blankdata %>%
  mutate(dimdv = ifelse(VenDor >= 13, 
                        "Dorsal", 
                        "Ventral"),
         dimml = ifelse(MedLat > 17, 
                        "Lateral",
                        ifelse(MedLat < 5,
                               "Medial",
                               NA)),
         dim = ifelse(is.na(dimml), dimdv, dimml)) %>%
  rowwise() %>%
  mutate(symdim = ifelse(MedLat == round(symline$mlvals[which(symline$apvals == AntPos)]),"SymLine", dim)) %>%
  ungroup()

#blankOB
plot_ly(dvml_blank, x = ~MedLat, y = ~AntPos, z = ~VenDor, color = ~VenDor, 
        text = ~paste('Dim: ', dim,
                      '<br>AP:', AntPos, 
                      '<br>ML:', MedLat,
                      '<br>VD:', VenDor),
        marker = list(size = 8, symbol = "circle"),
        type = 'scatter3d',
        mode = 'markers',
        hoverinfo = "none",
        hovertext="none") %>%
  layout(scene = list(xaxis = list(title = 'Medial-Lateral',
                                   backgroundcolor="rgb(200, 200, 230",
                                   gridcolor="rgb(255,255,255)",
                                   showbackground=TRUE,
                                   zerolinecolor="rgb(255,255,255"),
                      yaxis = list(title = 'Anterior-Posterior',
                                   backgroundcolor="rgb(200, 200, 230",
                                   gridcolor="rgb(255,255,255)",
                                   showbackground=TRUE,
                                   zerolinecolor="rgb(255,255,255)"),
                      zaxis = list(title = 'Ventral-Dorsal',
                                   backgroundcolor="rgb(200, 200, 230",
                                   gridcolor="rgb(255,255,255)",
                                   showbackground=TRUE,
                                   zerolinecolor="rgb(255,255,255)")))
#rankedolfr2 is a tidy_result ranked tibble using all intact olfr tpm
#rankedchemo is a tidy_result ranked tibble using all intact chemosensory gene tpm
#orbyorcor <- tibble(olfrname = ranked_olfr2$olfrname, orp50 = ranked_olfr2$p50, chemp50 = ranked_chemORs$p50) %>% group_by(olfrname) %>% summarise(cor = cor(orp50, chemp50)) %>% left_join(info, by = "olfrname")
orXor_cor <- read_csv("~/Desktop/obmap/r_analysis/3dimOB/output/chemoVolfr_p50.csv")
## Parsed with column specification:
## cols(
##   olfrname = col_character(),
##   cor = col_double(),
##   class = col_double(),
##   tan_zone = col_character(),
##   tzsimple = col_double(),
##   oe_region = col_character(),
##   RTP = col_character(),
##   known = col_logical(),
##   lowTPM = col_logical()
## )
ggplot(orXor_cor) + geom_violin(aes(oe_region, cor)) + 
  ggtitle("150iter, using all Chemo does not significantly alter voxel probability rankings")

#load tidy result and ranked

#load model output
tidy_result <- readRDS("~/Desktop/obmap/r_analysis/3dimOB/tidyresults/200i_allchemo_dimrep12_allvox.RDS")

ranked <- tidy_result %>%
  arrange(gene) %>%
  mutate(olfrname = rep(colnames(kzmY), each=max(tidy_result$voxel))) %>%
  group_by(olfrname) %>%
  mutate(voxrankperOR = rank(desc(p50))) %>%
  ungroup() %>%
  group_by(voxel) %>%
  mutate(ORrankpervox = rank(desc(p50))) %>%
  ungroup() %>% group_by(gene) %>% 
  mutate(allVoxAvgp50 = mean(p50)) %>% 
  ungroup() %>% group_by(gene, AntPos) %>% 
  mutate(AP_Avgp50 = mean(p50), AP_SNR = p50/AP_Avgp50) %>%
  ungroup() %>% group_by(gene, MedLat) %>% 
  mutate(ML_Avgp50 = mean(p50), ML_SNR = p50/ML_Avgp50) %>% 
  ungroup() %>% group_by(gene, VenDor) %>% 
  mutate(VD_Avgp50 = mean(p50), VD_SNR = p50/VD_Avgp50) %>% 
  ungroup() %>% 
  mutate(voxSNRdim = (AP_SNR + ML_SNR + VD_SNR)/3) %>% 
  group_by(gene) %>% 
  mutate(geneRankSNR = rank(desc(voxSNRdim))) %>% 
  ungroup() %>% group_by(voxel) %>% 
  mutate(voxRankSNR = rank(desc(voxSNRdim))) %>% 
  ungroup() %>% 
  select(p2.5, p50, p97.5, AntPos:ORrankpervox, geneRankSNR, 
         voxRankSNR, AP_Avgp50:voxSNRdim, voxel)

3D plots of OR probability as selected and constrained using my own functions

Plot X best voxels for an OR as ranked by SNR

Each voxel holds a probability for every OR, as determined by the expression of that ORs in the dimensional slices that intersect at that voxel. The total of the probability for all ORs in a voxel sums to 1. SNR is signal-to-noise ratio, essentially how much does an individual voxel stand out from the average probability of its 3 dimensional slices (all voxels at that single dimension position). Idea was to eliminate ectopics whose high expression across the whole OB dominated the voxel probability for nearly all voxels. 72 is ~2.5% of total voxel positions and presents a robust visualization, meaning that increasing this number typically results in adding neighboring voxels.

Plot best X voxels based on raw p50 (probability)

Using raw probability here, its pretty much the same as SNR after removing ectopics.

Clustering top ranked points

Given a number of top probability ranked points, cluster spatial neighbors.

Cluster("Olfr881", 1000, minClustSize = 2, chooseOut = "plot", title = "Olfr881 - Clustered top 100 probability voxels")
Cluster("Olfr1377", 125, minClustSize = 2, chooseOut = "plot", title = "Olfr1377 - Clustered top 100 probability voxels")
# Cluster("Olfr16", chooseOut = 100, "plot")
# Cluster("Olfr17", chooseOut = 100, "plot")
# Cluster("Olfr15", chooseOut = 100, "plot")
# Cluster("Olfr160", chooseOut = 100, "plot")
# Cluster("Olfr155", chooseOut = 100, "plot")
# Cluster("Olfr1507", chooseOut = 100, "plot")

Pick the best Medial and Lateral cluster

Plot best Medial and best Lateral cluster from the set of clusters across a symmetry line that is calculated using single dimension heatmap data For clarity, it is not picking the 2 best clusters but the best cluster from the medial half and the best cluster from the lateral half.

BestML("Olfr881", clustersPerHalfBulb = 1, title = "Olfr881 - Highest Probability Medial and Lateral Clusters")
BestML("Olfr1377", clustersPerHalfBulb = 1, title = "Olfr1377 - Highest Probability Medial and Lateral Cluster")
# BestML("Olfr16", clustersPerHalfBulb = 1)
# BestML("Olfr17", clustersPerHalfBulb = 1)
# BestML("Olfr15", clustersPerHalfBulb = 1)
# BestML("Olfr160", clustersPerHalfBulb = 1)
# BestML("Olfr155", clustersPerHalfBulb = 1)
# BestML("Olfr1507", clustersPerHalfBulb = 1)

Check if new Dorsal constraint function improves output

New function checks if OR has evidence of dorsal OE expression (class 1, miyamichi zone < 2, matsunami diffE = dorsal). If so, check if either the Medial or Lateral predicted position is dorsal. If OR is likely dorsal and either Medial or Lateral glomerulus prediction is dorsal but other halfbulb glom is not dorsal, find a dorsal glom for that halfbulb In independent images for the more lateral glomerulus, 1377 seems slightly more anterior lateral than 881

wachtg <- c("Olfr881", "Olfr1377")
tgout <- ListDorML(wachtg, chooseOut = "data")
## [1] "Olfr881"
## [1] "NA 100"
## [1] "NA 300"
## [1] 500
## [1] 700
## [1] "Olfr1377"
## [1] 100
## [1] "Olfr881"
## [1] "NA 100"
## [1] "NA 200"
## [1] "NA 300"
## [1] "NA 400"
## [1] 500
## [1] 600
tgout
## # A tibble: 72 x 39
##       p2.5     p50   p97.5 AntPos MedLat VenDor olfrname voxrankperOR
##      <dbl>   <dbl>   <dbl>  <dbl>  <dbl>  <dbl> <chr>           <dbl>
##  1 0.00645 0.00665 0.00695     14      1     15 Olfr881             6
##  2 0.00583 0.00611 0.00644     11      1     15 Olfr881            11
##  3 0.00515 0.00551 0.00592     12      1     15 Olfr881            14
##  4 0.00476 0.00507 0.00544     14      1     16 Olfr881            18
##  5 0.00428 0.00474 0.00522     14      1     14 Olfr881            26
##  6 0.00428 0.00461 0.00491      9      1     16 Olfr881            29
##  7 0.00425 0.00458 0.00489     11      1     16 Olfr881            32
##  8 0.00408 0.00435 0.00462     11      1     17 Olfr881            40
##  9 0.00411 0.00434 0.00465      9      1     17 Olfr881            41
## 10 0.00387 0.00418 0.00457     12      1     16 Olfr881            46
## # … with 62 more rows, and 31 more variables: ORrankpervox <dbl>,
## #   rankofrank <int>, rawclust <int>, clustmaxp <dbl>, clustminp <dbl>,
## #   clustmeanp <dbl>, n <int>, clustmaxprank <int>, clustmeanprank <int>,
## #   clustsizerank <int>, clust_unique <int>, isCluster <dbl>, meanML <dbl>,
## #   meanAP <dbl>, meanVD <dbl>, side <chr>, class <dbl>, tan_zone <chr>,
## #   tzsimple <dbl>, oe_region <chr>, RTP <chr>, known <lgl>, lowTPM <lgl>,
## #   dorsalRating <dbl>, MedRank <int>, LatRank <int>, TopStep <dbl>,
## #   minSize <dbl>, clustIn <dbl>, sideRank <int>, test <chr>
tgplot <- tgout %>% filter(p50 == clustmaxp) %>% 
  select(AntPos, MedLat, olfrname) %>% 
  filter(AntPos > 10) %>% filter(MedLat > 10)
blankout <- dvml_blank %>% select(AntPos:MedLat) %>% 
  mutate(olfrname = "aaa") %>% unique()
doubleout <- bind_rows(tgplot, blankout)

ggplot(blankout) + geom_point(aes(AntPos, MedLat, color = olfrname), size = 3) +
  geom_point(data = tgplot, aes(AntPos, MedLat, color = olfrname), size = 5) +
  theme_cowplot() +
  theme(legend.position = "none")

#DorsalML(wachtg)

Find positions for 50 ORs enriched in Wachowiak Functional Imaging surface samples

Shawn provided samples from 8 OBs from 4 mice, each OB cut into 2 pieces with 1 piece representing the functional imaging surface. Using the above algorithm for picking the best medial and lateral cluster for a given OR. Need to update list for newest alignment updates and to include all FI ORs now that code was improved in terms of speed.

olfr_result <- read_csv("~/Desktop/rproj/obmap_inactive/wach_diffe/starrsem_aligned/out/wach_v16model_top25FIenriched.csv") %>% mutate(Used1 = ifelse(Used == 1, 1, 0))

func_sig <- olfr_result %>% filter(logFC > 0) %>% filter(FDR < 0.05) %>% mutate(pseudo = str_detect(Gene_name, "-ps")) %>% filter(pseudo == F) %>% arrange(FDR)
func_sig_olfr <- func_sig$Gene_name
func_sig_olfr
##  [1] "Olfr376"  "Olfr629"  "Olfr57"   "Olfr1046" "Olfr1122" "Olfr937" 
##  [7] "Olfr994"  "Olfr228"  "Olfr570"  "Olfr558"  "Olfr969"  "Olfr402" 
## [13] "Olfr683"  "Olfr597"  "Olfr19"   "Olfr578"  "Olfr550"  "Olfr566" 
## [19] "Olfr1496" "Olfr974"  "Olfr478"  "Olfr677"  "Olfr147"  "Olfr561" 
## [25] "Olfr957"  "Olfr971"  "Olfr633"  "Olfr231"  "Olfr1032" "Olfr922" 
## [31] "Olfr1023" "Olfr635"  "Olfr1031" "Olfr690"  "Olfr1010" "Olfr1019"
## [37] "Olfr609"  "Olfr874"  "Olfr510"  "Olfr1134" "Olfr160"  "Olfr197" 
## [43] "Olfr467"  "Olfr150"  "Olfr935"  "Olfr64"   "Olfr5"    "Olfr506" 
## [49] "Olfr1377" "Olfr338"  "Olfr1086" "Olfr1020" "Olfr202"  "Olfr1339"
## [55] "Olfr225"  "Olfr691"  "Olfr895"  "Olfr20"   "Olfr982"  "Olfr146" 
## [61] "Olfr1328" "Olfr1129" "Olfr432"  "Olfr152"  "Olfr1449" "Olfr490" 
## [67] "Olfr557"  "Olfr51"   "Olfr488"  "Olfr1448" "Olfr1154" "Olfr1128"
## [73] "Olfr1511" "Olfr881"  "Olfr133"
ggplot(olfr_result) + 
  geom_point(aes(logFC,-log10(FDR), alpha = 0.25, color = as.factor(Used1), size = 1.3)) +
  geom_vline(xintercept = 0) + 
  geom_hline(yintercept = -log10(0.05)) + 
  theme_cowplot() + 
  theme(legend.position = "none") + 
  xlab("nonFIsurface  <<<  log2FoldChange  >>>  FIsurface")

Plot 30 of the “best” FI surface ORs

Best in this case refers to highestFDR (aka how consistency enriched in functional imaging surface). Perhaps color by an adjusted FDR?

#plot only highest p50 voxel of each cluster
ListML(func_sig_olfr[1:30], chooseOut = "point", title = "Medial/Lateral cluster points for the top 30 FI surface enriched ORs")

#new dorsal constraint function
ListDorML(func_sig_olfr[1:30], chooseOut = "point", title = "DorFunc Medial/Lateral cluster points for the top 30 FI surface enriched ORs")

Class 1 vs Class 2 positions

Expect Class 1 OR positions to be primarily dorsal-anterior to dorsal-central. Filtered out some ORs based on max to mean TPM ratio lower than 10 (968 ORs remaining). These could be ORs that are poorly enriched/dropout or have non-traditional expression across the OB (requires closer investigation). A cutoff of 10 could be high but current runtime for 968 ORs is ~40mins, need to c++ the pairwise clustering or alter step based cluster finding.

#get summary statistics for each OR
kzmY[1:5,1:5]
##         Olfr299 Olfr109 Olfr281 Olfr1015 Olfr1347
## sample1    0.00    0.00    0.33     5.77      0.0
## sample2  277.35    0.00  602.80   228.19      0.3
## sample3    0.00   99.66 2451.54     0.11      0.0
## sample4    0.00    0.00 1882.25   340.79      0.0
## sample5    0.00    0.00 2612.72     0.00      0.0
max_tpm <- vector(mode = "numeric", length = ncol(kzmY))
mean_tpm <- vector(mode = "numeric", length = ncol(kzmY))
for (i in 1:ncol(kzmY)) {
  max_tpm[i] <- max(kzmY[,i])
  mean_tpm[i] <- mean(kzmY[,i])
}
ornames <- colnames(kzmY)

#Removing 114 ORs whose max was less than 10 times that of their mean across all samples
metrics <- tibble(ornames, max_tpm, mean_tpm) %>%
  rowwise() %>%
  mutate(max2mean = max_tpm/mean_tpm,
         lowmax2mean = ifelse(max2mean < 10, T, F))

goodORs <- metrics %>%
  filter(lowmax2mean == F) %>%
  select(ornames) %>%
  as_vector() 
length(goodORs) #972
## [1] 1242
#good_out <- ListDorML(goodORs, chooseOut = "data")
#write_csv(good_out, "~/Desktop/obmap/r_analysis/3dimOB/output/goodchemo_ldML_200925.csv")
good_out <- read_csv("~/Desktop/obmap/r_analysis/3dimOB/output/goodchemo_ldML_200925.csv")

#removed 4 ORs had NAs for class and oe_region
good_point <- good_out %>% 
  filter(p50 == clustmaxp) %>% 
  select(olfrname, AntPos, VenDor, MedLat, class, oe_region, tzsimple, side) %>%
  mutate(class_fct = as_factor(class)) %>%
  filter(!is.na(class)) %>%
  filter(!is.na(oe_region))

write_csv(good_point, "~/Desktop/obmap/r_analysis/3dimOB/output/goodpointORs972_allchemo_listdorML_200820.csv")

shell_point <- BestML("Olfr10", topMin = 100, topMax = 100, 
                    topBy = 0, chooseOut = "notbest") %>%
  select(olfrname, AntPos:VenDor) %>%
  mutate(class_fct = NA,
         oe_region = NA,
         tzsimple = NA)

#lets look at class 1 vs class 2 ORs, note that is it possible for a voxel to hold multiple OR cluster points
plot_ly(type = "scatter3d", mode = "markers") %>% 
  add_trace(data=shell_point, x=~AntPos, y=~MedLat, z=~VenDor, 
            color="shell", opacity=0.15,
            text = ~paste('Gene: ', olfrname, 
                          '<br>AntPos: ', AntPos, 
                          '<br>MedLat: ', MedLat,
                          '<br>VenDor: ', VenDor),
            marker = list(size = 6)) %>%
  add_trace(data=good_point, x=~AntPos, y=~MedLat, z=~VenDor, color=~class_fct,
            text = ~paste('Gene: ', olfrname,
                          '<br>Class: ', class,
                          '<br>AntPos: ', AntPos, 
                          '<br>MedLat: ', MedLat,
                          '<br>VenDor: ', VenDor),
            marker = list(size = 6, line = list(color = 'black', width = 0.5))) %>%
  layout(title = "Class(1/2) of cluster points for 968 ORs",
         scene = list(xaxis = list(title = 'Anterior-Posterior'),
                      yaxis = list(title = 'Medial-Lateral'),
                      zaxis = list(title = 'Ventral-Dorsal')))
#a look at proportions to deal with point density
class_props <- good_point %>% 
  group_by(olfrname) %>%
  mutate(VDmeanpos = round(mean(VenDor))) %>%
  ungroup() %>%
  mutate(isc1 = ifelse(class == 1, T, F)) %>%
  group_by(VDmeanpos) %>%
  summarise(count = n(),
            proportion_class1 = sum(isc1)/count)

class_props %>%
  ggplot() + 
  geom_bar(aes(VDmeanpos, proportion_class1), stat = "identity") + 
  ggtitle("Proportion of Class 1 ORs across Ventral-Dorsal Axis") + 
  xlab("Ventral  <<<    100um sections   >>>  Dorsal")

#which class1 ORs are in ventral sections (VD < 8)
c1ventral <- good_point %>%
  filter(class == 1) %>%
  filter(VenDor <= 8)

OE region positions (zonal expression of OR as determined by Matsunami Lab DiffE)

3 samples of dorsal OE vs 3 samples of ventral OE Could also examine relation to more discrete tan et al. zone indices but this is more readable.

#matsunami diffe oe
plot_ly(type = "scatter3d", mode = "markers") %>% 
  add_trace(data=shell_point, x=~AntPos, y=~MedLat, z=~VenDor, 
            color="shell", opacity=0.15,
            text = ~paste('Gene: ', olfrname, 
                          '<br>AntPos: ', AntPos, 
                          '<br>MedLat: ', MedLat,
                          '<br>VenDor: ', VenDor),
            marker = list(size = 6)) %>%
  add_trace(data=good_point, x=~AntPos, y=~MedLat, z=~VenDor, color=~oe_region,
            text = ~paste('Gene: ', olfrname,
                          '<br>OE Zone: ', oe_region,
                          '<br>AntPos: ', AntPos, 
                          '<br>MedLat: ', MedLat,
                          '<br>VenDor: ', VenDor),
            marker = list(size = 6, line = list(color = 'black', width = 0.5))) %>%
  layout(title = "OE Zone of cluster points for 968 ORs",
         scene = list(xaxis = list(title = 'Anterior-Posterior'),
                      yaxis = list(title = 'Medial-Lateral'),
                      zaxis = list(title = 'Ventral-Dorsal')))
#a look at proportions to deal with point density
good_point %>% 
  group_by(olfrname) %>%
  mutate(VDmeanpos = round(mean(VenDor))) %>%
  ungroup() %>%
  mutate(isdor = ifelse(oe_region == "Dorsal", T, F)) %>%
  group_by(VDmeanpos) %>%
  summarise(count = n(),
            proportion_dorsal = sum(isdor)/count) %>%
  ggplot() + 
  geom_bar(aes(VDmeanpos, proportion_dorsal), stat = "identity") + 
  ggtitle("Proportion of Dorsal OE Zone ORs across Ventral-Dorsal Axis") + 
  xlab("Ventral  <<<    100um sections   >>>  Dorsal")

tan_point <- good_point %>%
  mutate(tzsimplest = ifelse(tzsimple <= 5, 6- tzsimple, NA)) %>%
  filter(!is.na(tzsimplest))

#tanzone
plot_ly(type = "scatter3d", mode = "markers") %>% 
  add_trace(data=shell_point, x=~AntPos, y=~MedLat, z=~VenDor, 
            color="shell", opacity=0.15,
            text = ~paste('Gene: ', olfrname, 
                          '<br>AntPos: ', AntPos, 
                          '<br>MedLat: ', MedLat,
                          '<br>VenDor: ', VenDor),
            marker = list(size = 6)) %>%
  add_trace(data=tan_point, x=~AntPos, y=~MedLat, z=~VenDor, color=~tzsimplest,
            text = ~paste('Gene: ', olfrname,
                          '<br>Tan Zone: ', tzsimple,
                          '<br>AntPos: ', AntPos, 
                          '<br>MedLat: ', MedLat,
                          '<br>VenDor: ', VenDor),
            marker = list(size = 6, line = list(color = 'black', width = 0.5))) %>%
  layout(title = "Tan Zone of cluster points for 968 ORs",
         scene = list(xaxis = list(title = 'Anterior-Posterior'),
                      yaxis = list(title = 'Medial-Lateral'),
                      zaxis = list(title = 'Ventral-Dorsal')))
#Mayra/Antonio/Luis topics
topics <- read_csv("~/Desktop/obmap/r_analysis/3dimOB/input/ORsDegreesOfBelonging.csv") %>% rename(olfrname = X1)

topic_point <- good_point %>% 
  left_join(topics, by = "olfrname") %>%
  filter(!is.na(maxTopic)) %>%
  mutate(maxTopicFct = as_factor(maxTopic))

#1 is all over
#2 is dorsal, class 1 enriched (67 class 1 ORs, 99 class 2 ORs)
#3 is ventral
#4 is only 1 Olfr338, dorsal
#5 is posterior
plot_ly(type = "scatter3d", mode = "markers") %>% 
  add_trace(data=shell_point, x=~AntPos, y=~MedLat, z=~VenDor, 
            color="shell", opacity=0.15,
            text = ~paste('Gene: ', olfrname, 
                          '<br>AntPos: ', AntPos, 
                          '<br>MedLat: ', MedLat,
                          '<br>VenDor: ', VenDor),
            marker = list(size = 6)) %>%
  add_trace(data=topic_point, x=~AntPos, y=~MedLat, z=~VenDor, color=~maxTopicFct,
            text = ~paste('Gene: ', olfrname,
                          '<br>Tan Zone: ', tzsimple,
                          '<br>AntPos: ', AntPos, 
                          '<br>MedLat: ', MedLat,
                          '<br>VenDor: ', VenDor),
            marker = list(size = 6, line = list(color = 'black', width = 0.5))) %>%
  layout(title = "MaxTopic of cluster points for 583 ORs",
         scene = list(xaxis = list(title = 'Anterior-Posterior'),
                      yaxis = list(title = 'Medial-Lateral'),
                      zaxis = list(title = 'Ventral-Dorsal')))
#visualize topic 1 ORs, color indicates % belonging
plot_ly(type = "scatter3d", mode = "markers") %>% 
  add_trace(data=shell_point, x=~AntPos, y=~MedLat, z=~VenDor, 
            color="shell", opacity=0.15,
            text = ~paste('Gene: ', olfrname, 
                          '<br>AntPos: ', AntPos, 
                          '<br>MedLat: ', MedLat,
                          '<br>VenDor: ', VenDor),
            marker = list(size = 6)) %>%
  add_trace(data=topic_point %>% filter(maxTopicFct == 1), x=~AntPos, y=~MedLat, z=~VenDor, color=~T1,
            text = ~paste('Gene: ', olfrname,
                          '<br>Tan Zone: ', tzsimple,
                          '<br>AntPos: ', AntPos, 
                          '<br>MedLat: ', MedLat,
                          '<br>VenDor: ', VenDor),
            marker = list(size = 6, line = list(color = 'black', width = 0.5))) %>%
  layout(title = "% Belonging Topic1 of cluster points for 583 ORs",
         scene = list(xaxis = list(title = 'Anterior-Posterior'),
                      yaxis = list(title = 'Medial-Lateral'),
                      zaxis = list(title = 'Ventral-Dorsal')))
good_point %>% 
  group_by(olfrname) %>%
  mutate(VDmeanpos = round(mean(VenDor))) %>%
  ungroup() %>%
  mutate(isdor = ifelse(oe_region == "Dorsal", T, F)) %>%
  group_by(VDmeanpos) %>%
  summarise(count = n(),
            proportion_dorsal = sum(isdor)/count) %>%
  ggplot() + 
  geom_bar(aes(VDmeanpos, proportion_dorsal), stat = "identity") +
  coord_flip() +
  theme_cowplot()
## `summarise()` ungrouping output (override with `.groups` argument)

#tanzone
plot_ly(type = "scatter3d", mode = "markers") %>% 
  add_trace(data=shell_point, x=~AntPos, y=~MedLat, z=~VenDor, 
            color="shell", opacity=0.1,
            text = ~paste('Gene: ', olfrname, 
                          '<br>AntPos: ', AntPos, 
                          '<br>MedLat: ', MedLat,
                          '<br>VenDor: ', VenDor),
            marker = list(size = 6),
            hoverinfo = "none",
            hovertext = "none") %>%
  add_trace(data=tan_point, x=~AntPos, y=~MedLat, z=~VenDor, color=~tzsimplest,
            marker = list(size = 8, line = list(color = 'black', width = 0.25)),
            hovertext = "none",
            hoverinfo = "none") %>%
  layout(title = "Tan Zone of cluster points for 968 ORs",
         scene = list(xaxis = list(title = 'Anterior-Posterior'),
                      yaxis = list(title = 'Medial-Lateral'),
                      zaxis = list(title = 'Ventral-Dorsal')))

Plot the 28 ORs enriched by 1% methyl-acetophenone

Color indicates logFC (pS6 strength of activity) Enriched defined as FDR < 0.05 and logFC > 1 pS6IP done by Maira. Control files are average no odor control metafiles (kevin 2020-04-23 created from starrsem of kzlabref).

Would be happy to create 2D plots or plot additional odors that Matt has imaged and we have corresponding pS6 data for if Hiro can identify shared odorants in both datasets.

plot heatmap peaks and calc dist to DorML

heatmap_peaks <- read_csv("~/Desktop/obmap/r_analysis/3dimOB/input/heatmap_peaks.csv")
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   olfrname = col_character(),
##   side = col_character(),
##   tan_zone = col_character(),
##   oe_region = col_character(),
##   RTP = col_character(),
##   known = col_logical(),
##   lowTPM = col_logical()
## )
## See spec(...) for full column specifications.
Heat3D("Olfr881", dimrep = 1, dv = "oe", raw = F, chooseOut = "plot")
all_heat_in_good <- heatmap_peaks$olfrname[which(heatmap_peaks$olfrname %in% good_point$olfrname)]
#ahig_dist <- DistHeat3D(all_heat_in_good)
#write_csv(ahig_dist, "~/Desktop/rproj/obmap/allmice/v21_gen25/heatmapORs_distance_to3D.csv")
ahig_dist <- read_csv("~/Desktop/obmap/r_analysis/3dimOB/output/heatmapORs_distance_to3D.csv")
## Parsed with column specification:
## cols(
##   olfrname = col_character(),
##   side = col_character(),
##   distance = col_double(),
##   ap_ml_vd_3d = col_character(),
##   ap_ml_vd_heat = col_character()
## )
ahig_dist %>% ggplot(aes(side, distance)) + geom_violin()

#TODO annotate positions of mombaerts 5OR mouse develop function that quantifies distance between plotbestML and assigned positions of mombORs define coordinates functional imaging surface define coordinates of zone boundaries develop function that quantifies for each OR, how many plotbestML gloms are in surface